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Abstract 

Background: Current research in network reverse engineering for genetic or metabolic networks very often does 
not include a proper experimental and/or input design. In this paper we address this issue in more detail and 
suggest a method that includes an iterative design of experiments based, on the most recent data that become 
available. The presented approach allows a reliable reconstruction of the network and addresses an important issue, 
i.e., the analysis and the propagation of uncertainties as they exist in both the data and in our own knowledge. 
These two types of uncertainties have their immediate ramifications for the uncertainties in the parameter 
estimates and, hence, are taken into account from the very beginning of our experimental design. 

Findings: The method is demonstrated for two small networks that include a genetic network for mRNA synthesis 
and degradation and an oscillatory network describing a molecular network underlying adenosine 3-5' cyclic 
monophosphate (cAMP) as observed in populations of Dyctyostelium cells. In both cases a substantial reduction in 
parameter uncertainty was observed. Extension to larger scale networks is possible but needs a more rigorous 
parameter estimation algorithm that includes sparsity as a constraint in the optimization procedure. 

Conclusion: We conclude that a careful experiment design very often (but not always) pays off in terms of 
reliability in the inferred network topology. For large scale networks a better parameter estimation algorithm is 
required that includes sparsity as an additional constraint. These algorithms are available in the literature and can 
also be used in an adaptive optimal design setting as demonstrated in this paper. 



Findings 

Background 

Recent research in Systems Biology shows that the con- 
cept of 'network' turns out to be highly relevant at all 
levels of biological complexity [1,2]. In a modeling ap- 
proach based on networks, the components of the system 
under consideration are represented as nodes in a graph 
and the interactions as edges that connect pairs of nodes. 
The functioning of a network is determined by its dy- 
namics, i.e., the evolution in time of the nodes. The nodes 
usually represent the concentrations of some species and 
their time development is regulated via interactions with 
other species. At the cellular level, regulatory and meta- 
bolic networks are out spoken examples, but also signal- 
ing pathways make use of the network concept. Typical 
examples at the highest complexity levels are predator- 
prey models and ecological food webs which are in use 
since long. It is remarkable that not only the network 
concept forms a unifying element across all biological 
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aggregation levels, but also the mathematics needed to 
describe the dynamics of the network nodes show great 
similarities [3]. The central role of network modeling in 
biology implies that network inference is one of the major 
challenges of modern biology. The ultimate aim of net- 
work inference is to deduce the structure of the network 
as accurately as possible from data obtained in the ex- 
perimental practice. To obtain the necessary experimen- 
tal information, one usually perturbs the network in 
some way, hoping that the induced response is useful to 
be explored in some network inference approach. How- 
ever, if this is not done in a structured way, the harvest 
might be disappointing. In this paper, we consider the 
question how we could effectively and fast learn as much 
as possible of the network structure by designing the ex- 
perimental setup in an optimal way. It should be realized 
that the question we put here in a network inference con- 
text, is in its generality not new, since similar research is 
at the core of the mathematical subdiscipline called' sys- 
tems theory', e.g., [4], and more specifically, in the subdis- 
cipline' system identification [5]. The fascinating aspect 
is that insights obtained in the engineering practice may 
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also be of great relevance in a life sciences setting. In an 
engineering context this kind of research is usually re- 
ferred to as reverse engineering'. In this paper we also 
aim at connecting the still too much separated realms of 
scientists in biology and engineering. 

Before getting into details, it is important to remark 
that network inference may aim at different levels of ac- 
curacy [6,7]. The lowest level approach leads to a rela- 
tively rough impression of the network structure. Given 
the nodes, one uses the data to find out whether a coup- 
ling exists between any pair of nodes. In mathematical 
terms, the network is represented as a graph with undir- 
ected edges and the inter action is not made explicit, e.g., 
in terms of an equation. An undirected edge indicates 
that the dynamics of the two connected nodes are corre- 
lated due to some interaction, but the character of this 
interaction is not pointed out in detail. For example, the 
first node could affect the second one, but the interaction 
could also be the other way around. This case is espe- 
cially relevant if the data contain steady state levels of the 
nodes obtained after perturbing the system a consider- 
able number of times, e.g., via a knocking out procedure 
in gene networks. For this case a number of statistical in- 
ference methods have been developed that yield undir- 
ected graphs as outcome [8-11]. The next level of 
accuracy is to strive for a directed graph, in which the 
edges are arrows. The direction of an arrow then indi- 
cates a causal relationship [12]. The most sophisticated 
level is to deduce more detailed information on the char- 
acter of the interactions. One then tries to answer ques- 
tions such as: Is it a promoting or an inhibiting 
interaction (in regulatory networks), is it a linearly in- 
creasing, saturating, or decaying interaction (in metabolic 
networks), etcetera [13,14]. In this paper we aim at the 
third level of accuracy by including the estimation of the 
strengths of interactions. A common procedure to gener- 
ate data is to perturb the system under consideration in a 
more or less random way. The change in the observed 
data is then analyzed to infer information about its struc- 
ture. In practice this may lead to very poor results, since 
it is not assured that the chosen perturbation is really 
useful for inference purposes. In this paper we propose 
to follow a much more advanced approach, in which the 
perturbations are designed such that they contain the op- 
timal amount of information, given the available data. 
This approach also starts with a more or less randomly 
chosen perturbation, but we show that the information 
from such a first step can be effectively used to design a 
second perturbation that yields enhanced insight in the 
network structure. So, we maximize in advance the infor- 
mation content of the second perturbation, given the 
insight deduced from the first perturbation. This 
maximization is based on improving the condition num- 
ber of the so-called Fischer Information Matrix of the 



system. If required, this procedure could be repeated to 
refine the inference results further. 

Our research focuses on networks whose dynamics are 
described in terms of ordinary differential equations 
(ODE). We assume that it is possible to observe the net- 
work while it develops in time, so that the data consists of 
time series of some or all of the network components. 
The general ideas of the present approach are widely ap- 
plicable and not at all restricted to a special type of net- 
work or level of complexity. For illustrational purposes we 
will make use of a gene is taken from [15]. This system 
has been used earlier in several studies to explain and il- 
lustrate the principles of other inference methods in sig- 
naling and gene networks [16-18] and turns out to be 
very useful for this purpose. Next, we apply our method 
to the Laub-Loomis model [19], which describes oscilla- 
tions in excitable cells of Dictyostelium, and show that 
our new procedure is very effective in reconstructing the 
underlying network, just by slightly perturbing system and 
observing its recovery to its natural oscillatory behavior. 

The paper is organized as follows. In the next section 
we will first introduce two motivating examples in which 
reconstruction of a directed graph representing the 
interaction nodes of the network is the primary goal. 
Then we will further elaborate on the details of adaptive 
optimal input design [20] and explain the method. 
Finally, the methodology will be applied to the two mo- 
tivating examples explained earlier and some results will 
be presented and discussed. 

Methods 

Two motivating examples 

In this section we introduce two biological systems that 
will be used later on to illustrate how our adaptive opti- 
mal design approach works in practice. For both models 
characteristic parameter values are taken from the litera- 
ture. These values are used to generate artificial data 
and to mimic real practice these data are spoiled by add- 
ing noise. The present aim is to show how the para- 
meters can be efficiently estimated from the data in an 
adaptive manner. 

mRNA synthesis and degradation Let us consider the 
well-known Kholodenko case as a motivating example 
[15]. It consists of a small four gene network and the 
gene activity is reflected by the synthesis and degrad- 
ation of mRNA as follows: 
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where x(t) stands for differentiation of x(t) with respect 
to time [h], and the synthesis and degradation rates are 
given by the following non-linear relationships: 
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The values of the parameters in these expressions are 
presented in Table 1. 

For simplicity we start by assuming that each time an 
experiment starts them RNA transcription and degrad- 
ation rates are in equilibrium, i.e. there is a steady state 
for all four states in the model. This means x\{t) — 
x 2 (t) = x 3 (t) = # 4 (£) =0. Furthermore, we assume, just 
as in the original paper, that we can influence the max- 
imum transcription rates {V is) i=l,. . .A} thereby silencing 
the synthesis of mRNA at will. In other words, the max- 
imum transcription rates represent an input signal that 
we can manipulate (or perturb) in such a way that infor- 
mation regarding the structure of the network is 
revealed through the stimulus-response data. An inter- 
esting (and challenging) aspect is that, in practice, it is 
not known what the exact values of the input signals are. 
Put differently, the (constant) input stimuli (or excitation 
signals) are in this case a natural part of the 
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Parameter values for the RNA model. These values were taken from [15]. 



identification problem and these values are included in 
the parameter estimation problem. The ODE system in 
equations (1)-(12) can be written in the general form 



x — f(x, 0, u) 



(13) 



Here, x is the state vector x = (xyX^x^x^) 1 , u is the 
vector containing the inputs, so u = (Vi^V^V^V^) 1 ^ ,0 
is the vector containing the remaining parameters in 
Table 1, and f is the vector valued function defined in 
equations (5)-(12). 

If we linearise the above equations and evaluate the 
resulting equations at the steady state x ss a linear state 
space model is obtained that reads as 



d(5x) 
dt 



A Sx(t) + B 8u 



(14) 



where Sx(t) and Su are now the deviations from the 
steady state x ss and the corresponding (steady) input u ss 
that maintains this steady state. The matrix A in equa- 
tion (14) contains the so-called interaction coefficients 
of the network and these correspond to a directed graph 
where there is an interaction from node / to node i if the 
corresponding matrix element ^ 0. We note that, of 
course, this linearisation is a local property that depends 
on our choice of the linearisation point in state space 
where we derive the linear system (14). However, it is 
important to realize that, if nodes i and j are not con- 
nected, the matrix A has zeros at positions ij and ji, re- 
gardless of the linearisation point that is used! The 
question of optimal experimental design now comes 
down to choosing the entries of the vector Su in such a 
way that matrix pair [A,B] can be estimated in an opti- 
mal manner. 

Laub-Loomis model-the no-input case A well-known 
example of an oscillating network in systems biology is 
Laub-Loomis' model of the molecular network under- 
lying adenosineS'-S'-cyclic monophosphate (cAMP) as 
observed in populations of Dyctyostelium cells [19]. The 
model incorporates changes in the activities of cAMP 
and consists of a set of seven ordinary differential equa- 
tions (one for each state) that reads 
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where the parameter values for the constants {k b 
i=l,...,14/ are given in Table 2. For these parameter 
values the system shows oscillating behavior. 

In the first motivating example we used the maximum 
synthesis rates as system inputs. Here, we assume that 
this oscillatory system can be perturbed by kicking it out 
its steady limit cycle. Since this cycle is asymptotically 
stable, the system will converge back to it. It is just the 
way in which this return happens, that provides us with 
the additional (and necessary) information to estimate 
the network structure. The linearised version of the 
Laub-Loomis model reads 



d(5x) 
dt 



A Sx(t) 



(22) 



If we allow a Dirac pulse, we simulate a sudden change 
in the state vector values so that the system is 'kicked' 
from its steady oscillation and needs to recover from this 
shock', returning (after some time) to its initial oscilla- 
tory behavior. Put differently, in mathematical terms we 
allow 'inputs' to be of the form bS (t), where b is a vec- 
tor with dimensions equal to the state vector x(t), so 
that 



d(5x) 
dt 



A Sx(t) + b S(t) 



(23) 



If we assume the perturbation Su to be a Dirac delta 
function (or better, a (^-distribution), the question of op- 
timal experimental design now becomes one of finding 
the values of the entries in the vector b that yield an op- 
timal estimate A, i.e., an estimate A with lowest possible 
uncertainty bounds. 

Parameter estimation- towards adaptive experimental 
design 

For a given choice of our input perturbation, e.g., Su in 
our first motivating example that contains four 
silencing-percentages for each of them RNA synthesis 
processes, we now perform a simple step response ex- 
periment for each of the four inputs {u t (t), i=l>. . .A} 
and this yields an observation record of mRNA concen- 
trations. In Figure 1 we see an example of such data that 
includes 20% error in the observed mRNA concentra- 
tions after a random relative perturbation (between-50% 
and 50%) of all input parameters {V is) i = 1,. . .A} in the 
original (non-linear) model. 

Table 2 Parameter values Laub-Loomis 
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Figure 1 The adaptive input design loop. Figure 2 represents the 
adaptive input design methodology as presented in this paper in a 
flow-chart. 



Parameter values for the cAMP model. These values were taken from [19]. 



Since there are four parameters {V is , i = 1,...A}, there 
are four 'inputs' at our disposal that we can use as per- 
turbation candidates. In the following we assume that 
four separate perturbations of all inputs are generated 
with random values between _50% and 50% of the steady 
state input values contained in the vector x ss (that corre- 
sponds to the initial steady state u ss ). 

Schmidt et al [16], inspired by the methods and results 
of [15], applied linear regression to data obtained by vir- 
tual experiments using the original system model as a 
simulation tool. Based on calculated deviations from the 
steady state x ss a linear regression formula was obtained 
with the coefficients of matrices A and B as unknowns. 
Schmidt successfully computed the interaction coeffi- 
cients, including the unknown input perturbations, from 
step-response data. For clarity it is noted that the step- 
response data represents the dynamic transient behavior 
from one steady state to another steady state. This is 
different from the approach presented by Steinke and 
co-workers [21], where an interesting statistical method- 
ology was applied with the same idea of iterative experi- 
mental design in mind, but for the steady states only, i.e. 
with no dynamic transient behavior included. 

It should further be noted that the regression method 
applied here only works for relatively small scale pro- 
blems in which the number of states is not too large. 
This is mainly because of the limitations in the applied 
parameter estimation method (see below for details) 
which easily introduces identifiability problems for large 
scale networks where the number of parameters to be 
estimated in the Jacobi matrix is simply too large. The 
disadvantage of ordinary least squares estimation, how- 
ever, can be remedied for larger scale networks by more 
advanced parameter estimation algorithms that include 
sparsity of the network as an additional constraint (see 
e.g. [22] for a convex linear-programming approach.) 
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A discretized version of the linear model (14) with 
sampling interval At=t/ <+ i_ t/< reads 



x(k + l) = 0x(k)+ru(k) 
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In summary the parameter estimation method pro- 
ceeds as follows: 

Step 1. Collect N observations of the state vector x in 
the original non-linear model (13) with a sampling inter- 
val At Subtract the steady state x ss from these 
observations. 

Step 2. Calculate the difference between consecutive 
time instances of the obtained measurements in step (1) 
and stack these differences in a (n x (N _ 1)) matrix M 
with n the dimension of the state vector x and N the 
number of data points. 

Step 3. Take M n as the 2 nd until the last column of M 
and take M G as the 1 st until the last but one column of 
M. Here, the subscripts n and o stand for new and old to 
denote a cause-effect relation between the matrices M n 
and M G . 

Step 4. Augment the matrix M G with a (1 x (N _ 1)) 
row of Is. This is to estimate the input sequence 



/Au(/<) = Tu(k),k = 0,. . .,N- I}. Denote the resulting aug- 
mented matrix with M a . 

Step 5. An estimate of the matrices O and T is now 
calculated via linear regression (see [16]) as 



(O f) =M n M/(M a M a 



(28) 



Step 6. Translate the discrete time result [0 f] back 
to continuous -time (see e.g. the standard textbook [23]), 
yielding the estimates A and B. 

At this point it should be mentioned that the familiar 
parameter estimation method in the above does not ex- 
plicitly include any analysis of the uncertainty in the par- 
ameter estimates obtained, although this is of course an 
important measure of the quality of the parameter esti- 
mates that have been obtained using the simple regres- 
sion formulae (28). Our goal in the current paper is to 
perform model based experiments that minimize the un- 
certainty bounds on the parameter estimates on-the-fly. 
To explain best what we mean with this the reader is re- 
ferred to Figure 2 in which the adaptive optimal input 
design methodology is depicted as an iterative loop. 
Starting with an initial random perturbation/experiment 
one moves on to a first parameter estimate [A, B] of the 
Jacobi matrices [A, B] following the 6-step procedure 
out lined in the above. Since at this point an approxi- 
mate model is available, one can progress by designing a 
new perturbation/experiment in such a way that a norm 
of the so-called Fisher Information Matrix (FIM) is opti- 
mized. The FIM incorporates the parametric output sen- 
sitivities JqZ Q \ , where y is the output vector which, in 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

time [hr] 

Figure 2 Time series data of a step response for the Kholodenko network model. Time series for the variables x v . . .,x 4 in equations (1)-(4) 
after a random perturbation of the parameter 1/1 s. We can clearly observe step-response curves in these data. 
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the two case studies discussed in this paper equals the 
state vector x since all states are assumed measurable 
with measurement error e(k), (k=l,. . .,iV), that is 
assumed to be sampled from a Gaussian white noise se- 
quence with co-variance matrix Q. We note that the 
parametric output sensitivity J^. e ^ can be calculated via 
the original state space model (13) (see e.g. [20]) through 
simple differentiation of the state equations with respect 
to the parameter vector 0. The FIM represents a meas- 
ure of the information content of the given input-output 
data with regard to the parameter values of the para- 
meters 6 in the model. We choose to minimize the so- 
called modified E-norm of the FIM [24], which is 
defined as the ratio of the maximum eigenvalue of the 
FIM by its minimum eigenvalue, i.e., 



F[t f -e 



dyf 



dyf 



d6(T-6)^ dO(T;6) 

max{Ap } 
Emod min{Af} 



di 



(29) 



(30) 



where tf marks the duration (final time) of our experi- 
ment and where we have denoted with 6 the fact that 
the FIM can only be calculated based upon the best 
knowledge we have available of the parameter values 
and not based upon the true values 0*. In practice this 
optimization can be performed crudely by evaluating a 
large number, say N = 1000, of random perturbations 
(through simulation of the approximate model (13) with 
the pair [A, B]) and evaluate the modified E-criterium 
for each of these perturbations. Once the N perfor- 
mances have been calculated we choose the perturbation 
that is associated with the minimum of the N modified 
E-criteria. For clarity we note that optimizing the FIM 
with regard to some norm can be seen as away to better 
condition the matrix inversion in equation (28). For in- 
formation rich experiments the matrix M a can be 
shown to be better conditioned and, hence, the inversion 
is less prone to numerical errors that can easily result in 
case of a badly condition matrix. 

Having thus progressed to the bottom-right circle in 
Figure 2, we continue to apply the 'optimal' perturbation 
to the real system, yielding an additional set of observa- 
tions that can now be used for a second parameter esti- 
mation exercise. We have now turned full circle and can 
start again with an improved estimate of the system 
[A, B], i.e. [A ,B ]. Several iterations of the loop in 
Figure 2 will now yield parameter estimates that have 
been inferred from an increasingly rich set of observa- 
tions. In other words, we have maximized the information 
content of the model parameters to the best of our 
(current) knowledge of the network as to converge rapidly 
to the true parameter values contained in the matrix-pair 



[A, B]. Through minimization of the uncertainty bounds 
on the parameter estimates via optimal input design a 
better conditioned parameter estimation problem is 
gained and this is highly preferred above the ad-hoc ap- 
proach where random perturbations are used as an ex- 
perimental 'design'. Although the above methodology 
(and more specifically the parameter estimates that are 
obtained using simple linear regression as in [16]) can 
only be used for relatively small networks, the adaptive 
input design approach can easily be extended to larger 
scale problems if other parameter estimation algorithms 
are utilized that can better handle larger scale networks 
(see e.g. [22] for a linear-programming approach to the 
parameter estimation problem for large networks). The 
message we would like to convey in the current paper is 
not about a parameter estimation algorithm, but about 
utilizing such an algorithm in the best possible way as to 
obtain reliable parameter estimates. 

Results 

Kholodenko case study 

On the basis of a random perturbation (with maximum 
perturbations of 50% for all four maximum transcription 
rates in the original non-linear model) we obtained arti- 
ficial time series for all four components of x, i.e. the 
four mRNA concentrations, with a time-interval AT = 1 
minute and for a length of 0.5hr, so that N = 30. We 
added 5% noise on these readouts. Furthermore, we 
assumed the matrix B to be diagonal. This implies that 
the perturbations Su influence the corresponding mRNA 
state directly. The first data sets were generated by per- 
forming 5 random perturbation experiments in a row 
and estimating the matrix A from these 5 experiments, 
yielding an estimate for A that is summarized in equation 
(33). The numerical values for A (deduced from 5 random 
experiments) together with their true values are: 



/-6.45 -2.92 
0.00 -8.17 
-2.31 2.80 
\ 0.00 0.00 



0.00 2.54 \ 
0.00 3.93 
-14.46 0.00 
10.22 -9.74/ 



/-7.15 -2.53 -1.69 3.20 \ 

-0.78 -8.19 0.039 5.83 

-1.51 1.63 -8.59 0.23 

V-1.74 1.65 0.37 -8.56/ 



(31) 



(32) 



The idea is now to try to improve the parameter esti- 
mates by performing an input design based on the 
model (13) with the initial estimates A and B that we 
obtained using the first data set-meaning that we used 
the first of the 5 random experiments to find a first esti- 
mate A for an input design for the second experiment. 
We then applied this procedure sequentially, i.e. the "op- 
timal" input found from the input design was fed to the 
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original non-linear system model to create a new data 
set and this data set was again used to estimate A, 
resulting in a second estimate of the Jacobi matrix, etc. 
This procedure was repeated until 5 experiments (using 
OID for the second to fifth experiment) were computed. 
The final values for the matrix entries in A were found 
to be: 

/-6.50 -3.14 1.36 2.43 \ 
2 _ 0.68 -7.66 3.20 3.15 ( . 

-3.05 3.91 -15.18 0.26 ^ ' 

\ 0.91 -1.45 10.25 -9.54 / 

These results have also been included in Figure 3 
which shows a comparison between the final estimate of 
the two approaches (OID and random). As to better get 
insight in the success rate of OID we repeated the above 
procedure 150 times. In Figure 4 a comparison of the 
two methods is presented in a graph to see the differ- 
ence in performance. Each point in this graph represents 
a pair of sum-of-squared-errors (of the parameters 
in the Jacobi matrix) for an optimal input design 
(x-coordinate) versus a random design (y-coordinate) 
series of 5 experiments in a row, repeated for 150 trials. 
The line y = x separates the two cases and we found that 
for 82% of the 150 pairs (OID, Random Design), the 
OID results were better, i.e. the x-coordinate was smaller 
than the y-coordinate. Apparently, random design can 
still be better than OID in 18% of the cases and this can 



be understood since a wrong first estimate A in the first 
OID experiment means that we design inputs based on 
erroneous information that points us in the wrong direc- 
tion for designing a better experiment. These results cor- 
roborate with the fact that if we reduce the 
measurement uncertainties on the sensors, OID 
becomes even better than Random Design because it has 
a better starting point in the series of 5 experiments. 
More specifically, after reducing the measurement error 
to only 1%, the success rate of OID versus Random was 
91%. 

Laub-Loomis case study 

The Laub-Loomis case study presents us with a far more 
challenging problem in comparison to the Kholondenko 
case study. This is mainly because in this particular case 
study 49 entries in the Jacobi matrix have to be esti- 
mated with linear regression, which is a large number. 
Furthermore, the non-linearity of the model introduces 
additional difficulties. It is therefore required to use rela- 
tively small perturbations in the initial condition x(0) as 
to not divert too far from the linearity assumption that 
underpins the linearized perturbation model. Since small 
perturbations are difficult to trace back if the noise on 
the data is too high, we are certainly limited in our 
possibilities. 

To perturb the oscillatory Laub-Loomis model we 
apply a shock'-effect to each of the entries in the state 
vector x(£) at time t = 0 and observe how the system 



Random Design Estimates 
Optimised Design Estimates 
True Values 



Figure 3 Parameter estimates drawn from a random perturbation experiment and an optimized perturbation experiment. Parameter 
estimates (and true values) for a random and an optimized experiment that was designed on the basis of the first estimate of A. The sixteen 
parameters are all entries in the Jacobian matrix A. We observe a substantial improvement of the parameter estimates once model based design 
was performed. 
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Sum of Squared Errors OID 

Figure 4 Comparison plot of Random Design versus Optimal Input Design (OID). This graph shows the sum-of-squared-errors in the 
parameters (Jacobi matrix entries) for 150 runs of our algorithm for a series of 5 experiments. Each run is a point in this graph with coordinates 
(Random, OID) where Random stands for the sum-of-squared-errors for the Random run and OID stands for the sum-of-squared-errors for the 
OID run that started with the same initial estimate of the matrix A that was obtained in the first random experiment. 

V . J 



recovers to its normal oscillatory behavior. An example 
of data from a random initial perturbation are presented 
in Figure 5. For the sampling interval of this observation 
record we assumed At to be 0.1 hr. Furthermore, the data 
set includes 20% artificial noise on the deviated observa- 
tions, i.e. the observed value x obs minus its nominal value 



x nom where the nominal value refers to an unperturbed 
oscillatory behavior. It can be observed in Figure 5 that a 
recovery of the initial perturbation to a steady oscillatory 
behavior (solid lines in figure) is present. 

We found in the initial run that the estimated param- 
eter values are still far off from the true parameter 




0 5 10 15 20 25 

time [hr] 

Figure 5 Response of the Laub-Loom is model for a random initial perturbation. Oscillatory behavior after a random initial perturbation. 
The solid lines represent unperturbed behavior while the dotted lines represent oscillations due to a small perturbation in the initial condition. 
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values. In order to find better estimates these initial esti- 
mates are now used for the next experimental design. 
Hereto the model (22) was evaluated for 250 random 
perturbations and for each of those 250 (virtual) experi- 
ments a modified E-criterion for the FIM was calculated. 
Taking the minimal value of these E-criteria lead to an 
improved perturbation that was subsequently applied to 
the real system (15) -(21). The second data-set was then 
used in a second parameter estimation exercise and OID 
procedure was repeated another time so that 3 experi- 
ments in total were used. In Figure 6 we see the final es- 
timation results after 1 random and 2 OID experiments. 
These are compared with a complete random experi- 
ment series of 3 experiments in total. Again, we see an 
improved parameter estimation result- almost all param- 
eter estimates are closer to their true values in compari- 
son with the random parameter estimation result. 

The above procedure of simulating three experiments 
and apply random design versus OID was then repeated 
100 times as to better compare the OID approach to a 
completely random approach. In Figure 7 we see a com- 
parison graph for pairs (OID, Random) together with 
the line y = x that represents the line of equal perform- 
ance. From the 100 trials (of three experiments each) we 
found that 74% of the OID results were better than a 
completely random design. In addition the average sum- 
of-squared-errors of the estimates of the entries in the 
Jacobi matrix were 214 for OID and 745 for the random 
experimental design, meaning that, on average, in this 
particular case OID performs three times better in terms 



25 p 
20 - 



of average uncertainty on the parameter estimates. This 
is certainly significant, even more after realizing that the 
Jacobi matrix is time-varying and 49 estimates need to 
be recovered from the noised at a provided. We also 
tried other noise-levels and sampling frequencies and 
found that when the noise on the data increases to 50% 
or more, there is no distinction between the random de- 
sign and OID visible. Clearly, the signal to noise ratio is 
so low in these cases that the parameter estimation algo- 
rithm cannot recover more information from the data 
on the basis of OID. But in such cases where the noise 
levels are just too high, clearly, there is not much to be 
gained from any approach since the information content 
of the data set is so low. 

Discussion and conclusions 

Adaptive optimal experimental design is a natural ap- 
proach for the recovery of the connections in an inter- 
action network based on a limited number of 
experiments. Clearly, the methodology as presented in 
this paper yields promising results as was demonstrated 
in two case studies. The idea we have pursued is that 
adaptive or sequential input design closes the model 
identification loop of experimentation and subsequent 
calibration of the model, intelligently taking in to ac- 
count the most recent knowledge that is available. More 
specifically, this means that the most recent parameter 
estimate of the Jacobi matrix pair [A, B] are used for 
subsequent analysis of the uncertainty propagation in 
the network. An essential feature proposed in our 



Random Design Estimates 
Optimised Design Estimates 
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Figure 6 Parameter estimates drawn from an optimized perturbation experiment. Parameter estimates (and true values) for a random 

experiment and an optimized experiment that was designed on the basis of the first estimate of A. The parameters are all entries in the Jacobian 

matrix A. We again observe improvement of the parameter estimates once model based design was performed. 
\ . ) 
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Sum of Squared Errors OID 

Figure 7 Comparison plot of Random Design versus Optimal Input Design (OID). This graph shows the sum-of-squared-errors in the 
parameters (Jacobi matrix entries) for 100 runs of our algorithm for a series of 3 experiments. Each run is a point in this graph with coordinates 
(Random, OID) where Random stands for the sum-of-squared-errors for the Random run and OID stands for the sum-of-squared-errors for the 
OID run that started with the same initial estimate of the matrix A that was obtained in the first random experiment. 



methodology (and also demonstrated in two case stud- 
ies) is that model based experiments are now promin- 
ently included in the loop and here with the best 
knowledge that is available, i.e. the most recent estimate 
[ A , B ] of the parameter vector 6. This yields a more 
carefully designed parameter estimation problem that 
allows a more reliable network reconstruction. 

One could argue that the linearized system introduces 
limitations since, of course, the matrices A and B de- 
pend on the linearization point chosen. Especially if the 
Jacobi matrix structure changes significantly (in terms of 
zeros and non-zeros) this can introduce a problem. A 
remedy would then be to repeat the algorithm for differ- 
ent values of the state vector x. However, we emphasize 
that if there is no connection between node i and ; in 
the network then regardless of the linearization point 
there will always be a zero in the if 1 entry. If, therefore, 
a zero is present in more than one linearization point, 
the chances are increasingly high that there is indeed no 
connection between node i and ; in the network struc- 
ture. Our approach contributes to the field of systems 
biology where the recovery of networks based on 
stimulus -response data is a central topic that has already 
caught a lot of attention [6,7]. In systems engineering, 
model based experimentation has matured substantially 
but here the design is usually developed for recovery of 
a set of model parameters in the original non-linear 
model structure and not for network inference, see e.g. 
[20] where the input design is formulated in a recursive 



manner, meaning that the data and subsequent planning 
of a new experiment are processed/analysed on-line. Of 
course, the results presented in the current paper do not 
include large scale networks with several hundreds of 
differential equations. But, then, this is not our main 
message. Rather, we have shown that OID based on 
dynamic time series of simple perturbation experiments 
leads to a better overall performance of the parameter 
estimation algorithm. In addition, the idea of adaptive 
input design may very well be combined with parameter 
estimation algorithms that are especially geared for such 
large scale problems and we think this will almost cer- 
tainly improve the efficiency of subsequent experimenta- 
tion and calibration as a means to unravel the 
underlying structure in a network model of a biological 
system. 
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